function [intensityPole, intensityMag, intensityT, intensityQ, intensityM] = FC_RadiativePower(file, namePool, freq)
    %FC_RADIATIVEPOWER 通过CST导出的表面电流计算辐射能(主要是基于磊哥的文件写的)
    %@since: 2022-04-27 10:22:30
    %        2022-05-23 20:47:01 更新位置
    %        2022-05-24 10:20:03 可以读取有无标题的两种文件类型
    %@update: 输入更新内容
    %@type:function
    %    *参数说明*
    %       @param [file] 文件路径
    %       @param [namePool] 命名的顺序
    %       @param [freq] 频率点的取值
    %       @return [返回名] 依次返回各个极子的能量，结构体形式存储
    %    *使用说明*
    %       [IP,IM,IT,QE,QM] = FC_RADIATIVEPOWER('\SC\',namePool,freq)
    %
    % ———————————————— POWERED BY YINWEI ————————————————
    c = 3e8; % 光速[m/s]
    w = 2 * pi * freq; % 角频率
    %% 导入数据，保存成cell
    for n = 1:length(namePool)
        disp(n);
        str = [file, num2str(namePool(n)), '.txt'];
        importedDataName{n}=importdata(str);

        if class(importedDataName{n}) == 'double'
            % txt 文件没有标题时，读取的类型直接就是一个cell
            importedDataValue{n} = importedDataName{n};
        elseif class(importedDataName{n}) == 'struct'
            % txt 文件有标题时，读取的类型是放入cell内部的struct
            importedDataValue{n} = importedDataName{n}.data;
        end

        kx{n} = complex(importedDataValue{n}(:, 4), importedDataValue{n}(:, 5));
        ky{n} = complex(importedDataValue{n}(:, 6), importedDataValue{n}(:, 7));
        kz{n} = complex(importedDataValue{n}(:, 8), importedDataValue{n}(:, 9));
        x = importedDataValue{1}(:, 1) * 1e-6; %x坐标
        y = importedDataValue{1}(:, 2) * 1e-6; %y坐标
        z = importedDataValue{1}(:, 3) * 1e-6; %z坐标
        s = importedDataValue{1}(:, 10) * 1e-12; %导出的面积
        Kxx(:, n) = [kx{n}]; %[kx{n}(:, 1)]
        Kyy(:, n) = [ky{n}]; %[ky{n}(:, 1)]
        Kzz(:, n) = [kz{n}]; %[kz{n}(:, 1)]
    end

    clearvars kx ky kz importedDataName importedDataValue str namePool

    xx = repmat(x, 1, n); % %x坐标重复扩充n列,[m]
    yy = repmat(y, 1, n); % %y坐标重复扩充n列,[m]
    zz = repmat(z, 1, n); % %z坐标重复扩充n列,[m]
    S = repmat(s, 1, n); % %面积s重复扩充n列,[m^2]
    r = [x, y, z]; % %坐标写成位矢,[m]
    R = repmat(r, 1, n); % %位矢r重复扩充n组,[m]
    clearvars x y z

    % ------------------------------
    %% 电偶极矩 已确认
    % Px = 1 ./ (i .* w) .* dot(S, Kxx);
    % Py = 1 ./ (i .* w) .* dot(S, Kyy);
    % Pz = 1 ./ (i .* w) .* dot(S, Kzz);
    Px = 1 ./ (1i .* w) .* sum(S .* Kxx, 1);
    Py = 1 ./ (1i .* w) .* sum(S .* Kyy, 1);
    Pz = 1 ./ (1i .* w) .* sum(S .* Kzz, 1);

    coefficient_P = ((2 .* w.^4) / (3 .* c^3));
    % %电偶极矩辐射能量
    Ipx = coefficient_P .* (abs(Px)).^2;
    Ipy = coefficient_P .* (abs(Py)).^2;
    Ipz = coefficient_P .* (abs(Pz)).^2;
    Ip = sqrt(abs(Ipx).^2 + abs(Ipy).^2 + abs(Ipz).^2);
    clearvars Px Py Pz coefficient_P
    intensityPole = struct('Ipx', Ipx, 'Ipy', Ipy, 'Ipz', Ipz, 'Ip', Ip);
    % clearvars  Ipx Ipy Ipz

    % ------------------------------
    %% 磁偶极矩 已确认
    % mx = (1 ./ (2 * c)) .* dot((conj(Kzz) .* yy - conj(Kyy) .* zz), S); %没有错误，dot共轭抵消
    % my = (1 ./ (2 * c)) .* dot((conj(Kxx) .* zz - conj(Kzz) .* xx), S);
    % mz = (1 ./ (2 * c)) .* dot((conj(Kyy) .* xx - conj(Kxx) .* yy), S);

    mx = (1 ./ (2 * c)) .* sum(((Kzz) .* yy - (Kyy) .* zz) .* S, 1);
    my = (1 ./ (2 * c)) .* sum(((Kxx) .* zz - (Kzz) .* xx) .* S, 1);
    mz = (1 ./ (2 * c)) .* sum(((Kyy) .* xx - (Kxx) .* yy) .* S, 1);

    coefficient_M = ((2 .* w.^4) / (3 .* c^3));
    Imx = coefficient_M .* (abs(mx)).^2; % %磁偶极矩辐射能量
    Imy = coefficient_M .* (abs(my)).^2;
    Imz = coefficient_M .* (abs(mz)).^2;
    Im = sqrt(abs(Imx).^2 + abs(Imy).^2 + abs(Imz).^2);
    clearvars mx my mz coefficient_M
    intensityMag = struct('Imx', Imx, 'Imy', Imy, 'Imz', Imz, 'Im', Im);

    % ------------------------------
    %% 环形偶极矩 已确认
    % rr = sqrt(x .* x + y .* y + z .* z); % %位矢的模
    rr = sqrt(sum(r.^2, 2));
    rr = repmat(rr, 1, n); % %rr重复扩充n列,[m]

    % Tx = (1 ./ (10 * c)) .* dot((xx .* (xx .* Kxx + yy .* Kyy + zz .* Kzz) - 2 .* rr.^2 .* Kxx), S);
    % Ty = (1 ./ (10 * c)) .* dot((yy .* (xx .* Kxx + yy .* Kyy + zz .* Kzz) - 2 .* rr.^2 .* Kyy), S);
    % Tz = (1 ./ (10 * c)) .* dot((zz .* (xx .* Kxx + yy .* Kyy + zz .* Kzz) - 2 .* rr.^2 .* Kzz), S);
    dotRJ = (xx .* Kxx + yy .* Kyy + zz .* Kzz); % 位移和电流点积
    cT = (1 ./ (10 * c));
    Tx = cT .* sum((xx .* dotRJ - 2 .* rr.^2 .* Kxx) .* S, 1);
    Ty = cT .* sum((yy .* dotRJ - 2 .* rr.^2 .* Kyy) .* S, 1);
    Tz = cT .* sum((zz .* dotRJ - 2 .* rr.^2 .* Kzz) .* S, 1);
    clearvars dotRJ cT
    % Tx = (1 ./ (10 * c)) .* dot((xx .* (xx .* Kxx) - 2 .* rr.^2 .* Kxx), S);
    % Ty = (1 ./ (10 * c)) .* dot((yy .* (yy .* Kyy) - 2 .* rr.^2 .* Kyy), S);
    % Tz = (1 ./ (10 * c)) .* dot((zz .* (zz .* Kzz) - 2 .* rr.^2 .* Kzz), S);

    coefficient_T = (2 .* w.^6) ./ (3 .* c^5);
    ITx = coefficient_T .* (abs(Tx)).^2; % %环形偶极矩辐射能量
    ITy = coefficient_T .* (abs(Ty)).^2;
    ITz = coefficient_T .* (abs(Tz)).^2;
    IT = sqrt(abs(ITx).^2 + abs(ITy).^2 + abs(ITz).^2);
    clearvars rr Tx Ty Tz coefficient_T
    intensityT = struct('ITx', ITx, 'ITy', ITy, 'ITz', ITz, 'IT', IT);

    % ------------------------------
    %% 电四极矩 已确认
    dotRJ = (xx .* Kxx + yy .* Kyy + zz .* Kzz); % 位移和电流点积
    cQe = (1 ./ (1i .* w));
    Qxx = cQe .* sum(((xx .* Kxx) + (xx .* Kxx) - 2/3 .* dotRJ) .* S, 1);
    Qxy = cQe .* sum(((xx .* Kyy) + (yy .* Kxx) - 2/3 .* dotRJ) .* S, 1);
    Qxz = cQe .* sum(((xx .* Kzz) + (zz .* Kxx) - 2/3 .* dotRJ) .* S, 1);
    Qyx = cQe .* sum(((yy .* Kxx) + (xx .* Kyy) - 2/3 .* dotRJ) .* S, 1);
    Qyy = cQe .* sum(((yy .* Kyy) + (yy .* Kyy) - 2/3 .* dotRJ) .* S, 1);
    Qyz = cQe .* sum(((yy .* Kzz) + (zz .* Kyy) - 2/3 .* dotRJ) .* S, 1);
    Qzx = cQe .* sum(((zz .* Kxx) + (xx .* Kzz) - 2/3 .* dotRJ) .* S, 1);
    Qzy = cQe .* sum(((zz .* Kyy) + (yy .* Kzz) - 2/3 .* dotRJ) .* S, 1);
    Qzz = cQe .* sum(((zz .* Kzz) + (zz .* Kzz) - 2/3 .* dotRJ) .* S, 1);
    clearvars dotRJ cQe

    coefficient_Qe = ((w.^6) / (5 .* c^5));
    IQxx = coefficient_Qe .* Qxx.^2; % %电四极矩辐射能量
    IQxy = coefficient_Qe .* Qxy.^2;
    IQxz = coefficient_Qe .* Qxz.^2;
    IQyx = coefficient_Qe .* Qyx.^2;
    IQyy = coefficient_Qe .* Qyy.^2;
    IQyz = coefficient_Qe .* Qyz.^2;
    IQzx = coefficient_Qe .* Qzx.^2;
    IQzy = coefficient_Qe .* Qzy.^2;
    IQzz = coefficient_Qe .* Qzz.^2;
    clearvars Qxx Qxy Qxz Qyx Qyy Qyz Qzx Qzy Qzz
    IQ = sqrt(abs(IQxx).^2 + abs(IQxy).^2 + abs(IQxz).^2 + abs(IQyx).^2 + abs(IQyy).^2 + abs(IQyz).^2 + abs(IQzx).^2 + abs(IQzy).^2 + abs(IQzz).^2);
    intensityQ = struct('IQxx', IQxx, 'IQxy', IQxy, 'IQxz', IQxz, ...
        'IQyx', IQyx, 'IQyy', IQyy, 'IQyz', IQyz, ...
        'IQzx', IQzx, 'IQzy', IQzy, 'IQzz', IQzz, ...
        'IQ', IQ);

    % ------------------------------
    %% 磁四极矩
    % RKx = conj(Kzz) .* conj(yy) - conj(Kyy) .* conj(zz); % %Rx(x,y,z)叉乘K(xx,yy,zz),x项;
    % RKy = conj(Kxx) .* conj(zz) - conj(Kzz) .* conj(xx); % %Rx(x,y,z)叉乘K(xx,yy,zz),y项;
    % RKz = conj(Kyy) .* conj(xx) - conj(Kxx) .* conj(yy); % %Rx(x,y,z)叉乘K(xx,yy,zz),z项;

    RKx = (Kzz) .* (yy) - (Kyy) .* (zz); % %Rx(x,y,z)叉乘K(xx,yy,zz),x项;
    RKy = (Kxx) .* (zz) - (Kzz) .* (xx); % %Rx(x,y,z)叉乘K(xx,yy,zz),y项;
    RKz = (Kyy) .* (xx) - (Kxx) .* (yy); % %Rx(x,y,z)叉乘K(xx,yy,zz),z项;

    Mxx = 1 ./ (3 .* c) .* sum(((RKx .* xx) + (RKx .* xx)) .* S, 1);
    Mxy = 1 ./ (3 .* c) .* sum(((RKx .* yy) + (RKy .* xx)) .* S, 1);
    Mxz = 1 ./ (3 .* c) .* sum(((RKx .* zz) + (RKz .* xx)) .* S, 1);
    Myx = 1 ./ (3 .* c) .* sum(((RKy .* xx) + (RKx .* yy)) .* S, 1);
    Myy = 1 ./ (3 .* c) .* sum(((RKy .* yy) + (RKy .* yy)) .* S, 1);
    Myz = 1 ./ (3 .* c) .* sum(((RKy .* zz) + (RKz .* yy)) .* S, 1);
    Mzx = 1 ./ (3 .* c) .* sum(((RKz .* xx) + (RKx .* zz)) .* S, 1);
    Mzy = 1 ./ (3 .* c) .* sum(((RKz .* yy) + (RKy .* zz)) .* S, 1);
    Mzz = 1 ./ (3 .* c) .* sum(((RKz .* zz) + (RKz .* zz)) .* S, 1);
    clearvars RKx RKy RKz

    coefficient_Qm = ((w.^6) / (20 .* c^5));
    IMxx = coefficient_Qm .* Mxx.^2; % %磁四极矩辐射能量
    IMxy = coefficient_Qm .* Mxy.^2;
    IMxz = coefficient_Qm .* Mxz.^2;
    IMyx = coefficient_Qm .* Myx.^2;
    IMyy = coefficient_Qm .* Myy.^2;
    IMyz = coefficient_Qm .* Myz.^2;
    IMzx = coefficient_Qm .* Mzx.^2;
    IMzy = coefficient_Qm .* Mzy.^2;
    IMzz = coefficient_Qm .* Mzz.^2;
    clearvars Mxx Mxy Mxz Myz Myy Myz Mzx Mzy Mzz coefficient_Qm
    IM = sqrt(abs(IMxx).^2 + abs(IMxy).^2 + abs(IMxz).^2 + abs(IMyx).^2 + abs(IMyy).^2 + abs(IMyz).^2 + abs(IMzx).^2 + abs(IMzy).^2 + abs(IMzz).^2);
    intensityM = struct('IMxx', IMxx, 'IMxy', IMxy, 'IMxz', IMxz, ...
        'IMyx', IMyx, 'IMyy', IMyy, 'IMyz', IMyz, ...
        'IMzx', IMzx, 'IMzy', IMzy, 'IMzz', IMzz, ...
        'IM', IM);
end

%% 绘图 %%总量

% YW_Plot_Fast({f / (1e12), Ip, f / (1e12), Im, f / (1e12), IT, f / (1e12), IQ, f / (1e12), IM}, {'-o', '-d', '-v', '-s', '-^b'}, ["\fontsize{28}\fontname{Helvetica}\itP", '\fontsize{28}\fontname{Helvetica}\itM', '\fontsize{28}\fontname{Helvetica}\itT', '\fontsize{28}\fontname{Helvetica}\itQ_e', '\fontsize{28}\fontname{Helvetica}\itQ_m'], ["Frequency (THz)", 'Radiation   power(a.u.)'])
